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Chapter 1 
INTRODUCTION 


Modal analysis has emerged as a valuable tool in many phases 
of the engineering design process. Complex vibration and acoustic 
problems In new designs can often be remedied through use of the method. 
Moreover, the technique has been used to enhance the conceptual under- 
standing of structures by serving to verify analytical models. 

In this thesis, a new modal parameter estimation procedure Is 
presented. The technique Is applicable to linear, time-invariant 
systems and accommodates multiple input excitations. In order to 
provide a background for the derivation of the method, section (1.1) 
briefly describes some modal parameter extraction procedures currently 
In use. Section (1.2) elaborates upon key features Implemented in the 
new technique. 

1. 1 Extraction Methods 

As better experimental measurement procedures and computation- 
al resources have become available, a need has likewise developed for 
more accurate and complete modal surveys. Quite sophisticated modal 
parameter extraction algorithms have consequently evolved. These 
methods include time and frequency domain analyses, and differ In the 
number and types of excitation applied to the structure. 

One principal category of extraction techniques recently 
Implemented in modal analysis is based upon representation of the 
response of a system in terms of characteristic "complex exponentials." 


1 
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In a general sense, the complex exponential techniques may be considered 
to Include the Ibrahim Time Domain (ITD) method, the Least Squares 
Complex Exponential (LSCE) method, and the Polyreference method, [1] In 
the ITD method, the free decay response of a structure is assumed to 
take the form 

2N X t 

x (t) - E 6 e r (1.1) 

r-1 ' 


In Eq. (1.1), x(t) is the response at N locations on the structure, <j> 
are modal vectors, and X are the eigenvalues of the system. By collec- 
ting the free decay response for several different initial conditions 
and curve-fitting the collected data to the form in Eq. (1.1), estimates 
of the modal parameters can be achieved. 

In the LSCE method, experimentally derived frequency response 
functions are fast Fourier transformed to impulse responses in the time 
domain. Analytically, these impulse responses are assumed to be 
expressible as 


h ij (t) 


F[H (w) ] 


N 

E 

r-1 


[A ijr e 


X t 
r 


- ^r C 
+ A. e ] 
i»r 


( 1 . 2 ) 


where FT ] and ( ) denote Fourier transformation and complex conjuga- 
tion. The quantities h^ (t) , (w) and A^^ in Eq. (1.2) are the 

impulse response, frequency response function, and rth complex residue 
for the ith and jth locations on the structure. As in the ITD method, 
estimates of the modal parameters are achieved by curve-fitting experi- 
mental data to the form of the complex exponentials in Eq. (1.2). 
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Both the LSCE and ITD algorithms have proven capable of 
generating accurate approximations of modal parameters. Yet, as with 
any method, criticisms of the approaches have arisen. One potential 
problem with the ITD method stems from the fact that it utilizes a 
single input model. Single input refers in this case to the particular 
initial condition selected in data acquisition. Reference [13] main- 
tains that the ITD method cannot identify repeated modes. 

In reference to the LSCE technique, the transformation of 
frequency response functions back to the time domain may introduce 
bias errors. Some difficulty has also been noted in determining an 
estimate of the number of active modes in a particular frequency range. 
[5] 

Unfortunately, the ITD And LSCE methods share an additional 
common disadvantage. While they both supply a consistent set of natural 
frequencies, neither approach generates a consistent set of modal 
vectors. [15] In this context, consistent means a unique set of 
calculated parameters. The LSCE method, for example, can supply a 
different set of modal vectors for data acquired using different exciter 
locations in a single run. [5] 

Techniques that generate a consistent set of eigenvalues and 
eigenvectors have been developed. The Polyreference method, as 
developed in (20) , is one such approach. Two other procedures that 
calculate consistent parameters are the Direct System Parameter 
Identification (DSPI) method derived by Leuridan [15] and the 
Simultaneous Frequency Domain technique of Coppolino. [7] The methods 
are similar in that they both approximate dynamic system matrices. 
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The DSPI algorithm is a multiple input, frequency domain 
parameter extraction technique. The approach may be formulated by 
transforming the governing equations 

[M]x(t) + [c]i(t) + [K]x(t) - f(t) (1.3) 

into their frequency domain form shown in Eq. (1.4). 

-w^[M]x(w) + jw[C]x(w) + [K]x(w) - f(w) (1.4) 

In Eqs. (1.3) and (1.4), fM] , fC] and [K] are N x N mass, damping, and 
stiffness matrices, respectively. The N-vectors x and f represent the 
response and excitation at N measurement locations on the structure. 
Because of the DSPI method seeks to estimate directly the entries in 
[M] , [C] and [K] , Eq. (1.4) is expanded for discrete frequencies and 
rearranged into a form suitable for calculation. 

[T]u - [V] (1.5) 

The matrix [T] in Eq. (1.5) is comprised of experimentally collected 
response, while [V] contains force spectra. The vector u is composed of 
the unknown elements in [M] , fC] and [K]. By controlling the location 
that the entries in Eq. (1.4) are indexed into the matrices in Eq. 
(1.5), constraints such as bandwidth and symmetry may be imposed upon 
the calculated [M] , [C] and [K] matrices. Actual solution of Eq. (1.5) 
proceeds using a least squares technique. Once the entries in the mass, 
damping and stiffness matrices are available, a general eigenvalue 
problem can be assembled and solved for the modal parameters. 
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While Che DSPI method does provide a sec of consistent modal 
parameters, the Indexing scheme used In assembling Eq. (1.5) Is quite 
complicated. Even more Importantly, If the number of measured responses 
is large, the DSPI technique requires an Inordinate amount of computer 
space. The method makes no provision for reducing the problem size. 

The SFD technique Is a single Input method that assumes that the 
response of 1 "Independent" locations on a structure is governed by the 
set of differential equations 


[M 1 ]x 1 (t) + fC 1 ]x 1 (t) + [K 1 ]x 1 (t) o [D ± Jf (t) (1.6) 

where 

[M^] ■ i x i mass matrix 

[C ± ] » 1 * i damping matrix 

[K^] ■ i x i stiffness matrix 

[D^] ■ 1 x 1 force distribution matrix 

f(t) ■ a single excitation time history 


As is explained shortly, the number of independent coordinates, 1, may 
be considerably less than the total number of measurement locations N. 
When Eq. (1.6) is Fourier transformed, multiplied by TH^] and 
rearranged, the fundamental equation of the SFD method emerges 


^(w) + j^[C i ][K 1 ][D t ]j 


x i (w) 


x jL (w) 


0 


(1.7) 


-f (w) 
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A 

the new matrices [C^], 

[K^] and 

A 

'V 

are 

defined In Eqs. 

(1.8A), (1.8B) , 

and ( 1 . 8C) . 







A 

[c i> 

- 

[MJ' 


(1.8A) 


[K i ] 

= 

[M ± ]* 

‘‘[Kj] 

(1.8B) 


A 

tv 

- 

[M ± ]' 

‘‘fV 

(1.8C) 


If response spectra or frequency response functions are collected, Eq. 
(1.7) can be expanded for each discrete frequency and solved for the 

A A A 

matrices [C^], [K^] and [D^] In a least squares sense. Once these 
matrices are obtained, a consistent set of eigenvalues and eigenvectors 
can be calculated from the standard eigenvalue problem 




(1.9) 


If the number of responses, N, Is quite large, the solution of 
Eq. (1.7) for all measurements, 1 “ N, may be cumbersome. The SFD 
algorithm however, provides a means of reducing the problem size If 
necessary. In the method, It Is possible to choose the 1 Independent 
locations to be a reasonably sized subset of the total N measurement 
locations. The remaining d - N - 1 locations are termed "dependent" 
coordinates. Solution of Eq. (1.7) using only the Independent coordi- 
nates Is then feasible. Because 1 < N, the eigenvectors calculated In 
Eq. (1.9) for the reduced problem are of a smaller dimension than full 
system eigenvectors. The entries of the full eigenvectors corresponding 
to dependent coordinates are calculated In a least squares sense from 
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the measured response spectra, or frequency response functions, and the 

eigenvectors for the Independent coordinates. 

To better Identify closely spaced modes, a version of the SFD 

algorithm has been derived that utilizes multiple excitation locations. 

Blair [3] has shown that minor modifications can result In significant 

improvement in the method's ability to resolve closely spaced modes. 

Essentially, the single excitation history f(t) in Eq. (1.6) is replaced 

with an N -vector of excitation time histories f(t), where N is the 
P - P 

number of exciters. As a result, Eq. (1.7) becomes 


x ± (w) + 


[jcjtKjHDjj 


x(w) 
x(w) 
-f (w) 


( 1 . 10 ) 


where [D. ] is now i x N . The solution of Eq. (1.10) for the matrices 
i P 

A A 

[C ± ] and [K^], and ultimately the calculation of system eigenvalues and 
eigenvectors, proceeds identically as in the SFD method. 

Selecting the independent coordinates used in model reduction 
poses a source of difficulty in the SFD technique and Blair's version of 
the method. In both cases, user judgment is required. In addition, 
each method uses the normal equations solution procedure, prone to 
numerical instabilities for poorly conditioned data. As a final point, 
the theoretical development of the SFD algorithm and Blair's version of 
the approach assumes that the system damping is proportional, a very 
limiting restriction. 
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1.2 Features of the Newly-developed Method 

Considering the relative merits of the techniques considered 
in section fl.l), the approach derived in this thesis is based upon 
Blair's modification of the SFD algorithm. The current study extends 
the work of Blair by investigating automatic procedures to reduce the 
effective problem size. It has been a primary intent of the research to 
evaluate means of decreasing the user interaction required in generating 
a reduced-order model. Two reduction strategies are outlined in the 
thesis. 

Furthermore, derivation of the parameter extraction equations 
proceeds without assumptions regarding the form of damping expected in 
the physical system. In contrast to the methods of Blair and Coppolino, 
a theoretical framework is presented that can admit skew-symmetric or 
non-proportional damping matrices. 

Lastly, the technique developed herein utilizes more stable 
solution procedures than those of Blair and Coppolino. A frequency 
response function synthesis technique is also described to facilitate 
modal parameter verification. 

The presentation of the new method begins in Chapter 2 with 
the description of the theoretical basis and computational organization 
of the technique. A method of verifying the calculated modal parameters 
is reviewed in Chapter 3. Chapter 4 examines transformations used in 
reducing large. Impractical problems to a reasonable size. Least 
squares solution procedures and applications to component mode synthesis 
are topics elaborated upon in Chapters 5 and 6, respectively. Both 
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analytical and experimental example problems are discussed in Chapter 7. 
The final chapter of the thesis contains conclusions and recommendations 
for further work. 



Chapter 2 

DEVELOPMENT OF THE PARAMETER EXTRACTION METHOD 


In developing a modal parameter estimation method for linear, 
time- invariant dynamic systems, it is assumed that there exists a 
discrete analytical model. Identical to that in Eq. (1.3), that 
accurately represents the dynamics of the structure. 


[M]x(t) + [c]x(t) + [K]x(t) - f(t) 


( 2 . 1 ) 


As in Chapter 1, x(t), x(t) and x(t) are the acceleration, velocity and 

displacement time histories at N discrete locations on the structure. 

Similarly, [M] , [C] and [K] are of N * N order and commonly referred to 

• as the mass, damping, and stiffness matrices, respectively. The 

N-vector f(t) contains input excitation functions, one for each of the N 
«* 

degrees of freedom. 

When dealing with such general systems, it is often more 
convenient to analyze the set of equations in Eq. (2.1) when they are 
organized in a different manner. In particular, a state vector 
formulation is frequently utilized which introduces a 2N-vector of 
unknowns X(t) . 


X(t) 


x (t) 
x(t) 


( 2 . 2 ) 


A close inspection of the alternate system 


[A]X(t) + [BlX(t) 


F(t) 


(2.3) 
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where 


[A] 


[0] [Ml 
[M] fC] 


(2.4) 


[B] 


-[M] [0] 

roj [k] 


(2.5) 


F(t) 


■ 

0 

f(t) 


( 2 . 6 ) 


reveals that Eqs. (2.1) and (2.3) are equivalent. The matrices [A] and 
[B] are of 2N * 2N order, while F(t) is a 2N-vector. 


2. 1 Model Reduction 

It has been noted in Chapter 1 that some methods are Ill- 
equipped to handle the numerous measurements that can result from 
practical experimental analyses. An analytical model such as Eqs. (2.1) 
or (2.3) that is large enough to represent all these experimental 
degrees of freedom may be difficult, or even impossible, to generate due 
to limited computer resources. Thus, the need to be able to reduce the 
size of the dynamic model becomes apparent. 

A common method used in analytical dynamics to approximate a 
large system of differential equations, as in Eq. (2.1), by a smaller 
system is Ritz analysis. [2, 8] Essentially, the solution x(t) is 
approximated by a linear combination of M Ritz basis vectors, such 


that 
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M 

x(t) 3 z y,(t)t|, (2.7) 

i-1 1 - 1 

where M < N. Because the approximation of x(t) is restricted to the 
subspace spanned by the Ritz vectors ip^, the accuracy of Eq. (2.7) 
depends upon the Ritz vectors chosen. Reasonable accuracy is achieved 
only if the Ritz vectors span a significant portion of the "motion 
space" of the system in Eq. (2.1). 

The fact that Ritz analysis can indeed decrease the dimension of a 
system of equations is apparent when Eq. (2.7) is written in matrix form 

x(t) s [ip 1 i(, 2 ... tf> M ]y(t) * [^yCt) (2.8) 

In Eq. (2.8), is an N x M matrix, and the original N-vector x(t) is 

approximated using a smaller order M-vector y(t). 

•* 

Ritz analysis may also be likened to another transformation 
method encountered in time series analysis. [4] The approach seeks to 
replace an unwieldy, large set of unknowns by a smaller, statistically 
similar set. For example, the technique may be applied in signal 
processing when it is necessary to transmit several signals over a 
smaller set of channels. To compensate for the limited number of 
channels, only a subset, or a reduced combination, of the original set 
are transmitted. The full, original set is then approximated upon 
reception of the subset of signals. [4] 

This process of reducing, and later reconstructing, a set of 
time histories is embodied in Eqs. (2.9), (2.10A) and (2.10B). 

y(t) - 


U c ]x(t) 


(2.9) 
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x(t) 

= [* R ]y(t) 

(2.10A) 

x(t) 

■ f'l' R ]y(t) + e(t) 

(2.10B) 


In equation (2.9), the original response N-vector x(t) is "condensed" 
via an M * N transformation matrix [ij^] to the reduced M-vector y(t). 

The condensed vector is assumed more convenient for processing, such as 
computation or transmission, than the original vector x(t). When 
processing of y(t) is complete, the original N-vector x(t) may be 
estimated from the M-vector y(t) using the N x M matrix O r ] in Eq. 
(2.10A). Equation (2.10B) simply identifies the error incurred in Eq. 
(2.10A) as an N-vector e(t). The notation in Eq. (2.10B) is convenient 
in later discussion of transformation methods. 

Comparison of Eqs. (2.8) and (2.10A) reveals that the recon- 
structing transformation is identical in form to a Ritz analysis trans- 
formation. As in the Ritz analysis, the accuracy of Eq. (2.10A) depends 

upon the selection of [# ] and Ctp_ ] . This selection process is dis- 

U K 

cussed in Chapter A. 

The above transformations can now be used to produce a reduced- 
order model. Substituting Eq. (2.10) into Eq. (2.1), and premultiplying 
by [<!»£,] yields a smaller set of governing equations in the new variable 
y(t) . 

[* c HMW R Jy(t) + r# c ncm R ]Y(t) + [* c ][KW R ]y(t) + - [<l> c ]f(t) 


( 2 . 11 ) 
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Neglecting the error terms and introducing a more compact notation 
results in the reduced system 

fM*b(t) + [C*] y(t) + TK*Mt) s [* c ]f(t) (2.12A) 

with 

[M*] - fd» c ] [MJ [ip R ] (2.12B) 

rc*] » [ d» c l [C] (2.12C) 

[K*] » [^ c l [K] ‘ (2.12D) 

where [M*], [C*], and [K*] are M * M mass, damping and stiffness matri- 
ces. Because of the omission of error terms, the equality in Eq. 

(2.12A) is necessarily approximate. However, further equations derived 
from Eq. (2.12A) are expressed as strict equalities. Lastly, it is 
worth noting that if ■> [t|^] T , then [M*], [ C* ] , and [K*] are identi- 

cal to the matrices obtained in a standard Ritz analysis reduction. 

[ 2 , 8 ] 


2.2 Frequency Domain Model 

With the advent of digital Fourier analyzers, frequency-domain 
analysis in structural dynamics has flourished. [22] These instruments 
provide a quick means of generating the frequency spectrum of a time- 
domain signal and may also directly compute frequency response func- 
tions. Because experimental frequency domain data is commonly available 
and may be easily transformed to represent accelerations, velocities or 
displacements, a frequency-domain equivalent of the reduced model in Eq. 
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(2.12) is sought. To this end, it is now necessary to Fourier transform 
Eq. (2.12). 


[M*]y(w) + [C*]Y (v) + [K*]y(w) s [* c ]f(w) (2.13) 

With T ] denoting a Fourier transform, the previous notation r.ay be 
defined 


Y(v) - F[Y(t)] 

(2.14A) 

Y(w) - F[Y(t)J 

(2.14B) 

Y (w) - F[Y(t)] 

(2. 14C) 


In practice, experimental acceleration spectra are more prevalent than 
velocity or displacement spectra. Using the property of Fourier 
transforms [10], ff^(t)] - (jw) n f(w), Eq. (2.13) can be written in 

terms of acceleration spectra alone. 

[M*Mw) + [C*] T^Y(w) + [K*]-^-Y(w) - r<J> c ]f(w) (2.15) 


Normally, a physical system undergoing multi-shaker modal 
testing has only a few points of excitation and many nodes. In this 
case, the full N-vector of nodal excitations, f(t), can be represented 
by 


f(t) - [D]p (t) (2.16) 


If N is the number of actual exciters, with N << N, then [D] is an 
P P 

N x N force distribution matrix and p(t) is an N -vector of excitation 
P - P 
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functions. When Eq. (2.16) Is Fourier-transformed and substituted Into 
Eq. (2.15), a new system results, 

[M*]y(w) + [C*] -^Y(w) + [K*] y(v) - [D*]p(w) (2.17) 

where [D*] = is an M x N matrix. Premultiplying by [M*l * 

C p 

produces the desired reduced-order frequency-domain analytical model for 
the system described originally by Eq. (2.1). 

Y(v) + [C] Y (w) [K] -^Y(w) - [D]p(w) (2.18) 

J - -v* 

A 

In this expression, [D] is an M x N matrix. 


2.3 Organization of Computations 

While the transformation scheme outlined in Eqs. (2.9) and 
(2.10) is defined for continuous, analytic functions, the approach is 
equally applicable to discrete data. In particular, the contracting 
transformation [i^] can be used to generate a condensed set of discrete 
experimental data for the reduced-order model in Eq. (2.18). If [x(w^)] 
represents measured vectors of response spectra, 

Tx(w )] - [x (w . ) x(w ) ... x (w )] (2.19) 

' ' ' N NxN 

w w 

then a reduced experimental set of spectra may be obtained via a form of 
Eq. (2.9) expanded for discrete measurements. 


[y(w ± ) 3 - 

MxN 

w 


... y(v )] 

- XT 


O c ] [x(w t )] 

MxN NxN 

w 


[y(w l ) y( w 2 ) 


( 2 . 20 ) 
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A quick inspection of Eqs. (2.9) and (2.10) reveals that results similar 
to Eqs. (2.19) and (2.20) may be written for the derivatives of x(t). 
Thus, 


ty(w 1 )] - [U» c Hx(w 1 )1 

MxN MxN NxN 
w w 


( 2 . 21 ) 


To accommodate the condensed experimental acceleration spectra 
obtained in Eq. (2.21), Eq. (2.18) must be expanded and rearranged for 
N w discrete frequencies. 


[C] [K] [D] 
Mx(2M+tT ) 


[-^2 r ( w i ) ] 

""i 


[-y(w.)1 


MxN 


w 


[-p(w t )l J (2M+N )XN 

p w 


( 2 . 22 ) 


Four new matrices are newly defined in Eq. (2.22). 




jw Y 1 (W 1 ) **• jw Y 1 (W N 5 
J 1 J N w 

w 


35T VV ... t y 2 (w n ) 
J 1 J N w 


7^ W • 


iw *M^ W N ^ 
J N w 

w 


MxN 


w 


(2.23A) 
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Y^Ol 

" W 1 


2 Y l^ w p *•* 


"* 2 Y 2 ^ 1 ^ * * * 

w. 


- — Y^w ) 

w., Nw 
N 

W 

- — 2 Y 2 (w n ) 


N 


v 


w 


“T W 
W 1 


2 Y M^ W N * 

W N W 

w 


MxN 


w 


[r<’< 1 >J 


W 


•• Y^Wj, ) 
w 


Y 2 ^ W J^ ••• Y 2 ^ w jqr ^ 


W ••• Y m (w N 5 


MxN 


w 


r-p 


-pjfWj) ... -P^Wjj ) 

w 


-P 2 (w 1 ) ... -P 2 (w n ) 


-P N ( w x ) •• 


• -Pn (u n > 

p w 


N x N 

p w 


(2.23B) 


(2.23C) 


(2.23D) 



If is given, the only unknown quantities in Eq. (2.22) are the 

A A A 

entries in [C], [K] and [D], These matrices implicitly represent the 
dynamic properties of the reduced physical system. A standard least- 

A A A 

squares solution for [C], [R] and [Dj may be formulated by transposing 
and separating the real and imaginary parts of Eq. (2.22). 



(2.24) 


Equation (2.24) may now be solved for [C], [K] and [D] by any one of 
several methods, two of which are described in Chapter 5. 

A A A 

Because the matrices [K], [C] and [D] are real "ahd constant, ! 
the time-domain equivalent of Eq. (2.18) is expressed 


Y(t) ♦ (C]y(t) + [Kfr(t) - [D]p(t) (2.25) 


In a manner analogous to that done for the full system, the state vector 
formulation for the reduced system in Eq. (2.25) may be written 


[I] [0] 
[0] [I] 


’ Y(t) ' 


' [ci m ' 


’ Y(t) ' 


' [DJ p(t) ' 

** 

+ 




- 


Y(t) 


[-i] roi 


Y(t) 


0 


(2.26) 


or more concisely 
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where 


Tit) 


y(t) 

y(t) 


(2.2TB) 


Although [A] Is a 2M x 2M Identity matrix, the present notation is 
retained to simplify discussion of synthesis in Chapter 3. 

A standard eigenproblem may be formulated for the eigenvalues 
and right eigenvectors of Eq. (2.27). If the time dependence of T(t) is 
given by r(t) ■ Xe , then the eigenproblem corresponding to Eq. (2.27) 
can be written 


A X + [B]X - 0 (2.28) 

r_r .r 

Equation (2.28) introduces a notational convention maintained throughout 
this thesis. All eigenvectors and static responses in the derivations 
that follow may be distinguished from other vectors in that they are 
expressed without an argument denoting time or frequency dependence. 

A A 

From Eq. (2.28), 2M eigenvalues A and 2M reduced eigenvectors, X f of 
length 2M may be extracted. The eigenvalues may be used directly as 
approximations to the eigenvalues A f of the original system (2.3). 

The reduced eigenvectors in Eq. (2.28) may be employed to 
achieve an approximation to the eigenvectors of the full system in Eq. 
(2.3). Recalling both the form of X(t) in Eq. (2.2) and the 
transformation in Eq. (2.10), we have 

x(t) 
x (t ) 


X(t) 


(2.29) 
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x(t) s 


(2.30) 


It is immediately apparent that 


X(t) s 




[ tp R l to 1 

fK. R ] Y (t) 


[0 ] [# R ] 


r(t) 


(2.31) 


Equation (2.31) defines the state vector form of the reconstructing 
2N x 2M matrix [f R ] . 


x(t) 3 [ ¥r ] r (t) 


(2.32) 




t* R ] [0 ] 

[0 ] r*, r ] 


(2.33) 


Finally, the transformation in Eq. (2.32) may be used to approximately 
reconstruct the full system eigenvectors from the reduced system 
eigenvectors. 


x s r>? D ]x 

r T R r 


(2.34) 



Chapter 3 

FREQUENCY RESPONSE FUNCTION SYNTHESIS 

Often in experimental modal analysis, it is necessary to 
verify that the set of modal parameters obtained is both accurate and 
complete over the frequency range of interest. A common means of 
validation is to generate an analytical set of frequency response 
functions, or FRF's, from the estimated modal parameters. The 
analytical FRF's may then be visually compared to experimentally 
collected FRF's. In this manner, spurious modes may be rejected from a 
set of calculated modal parameters. 


3. 1 Preliminaries 

The derivation of frequency response synthesis formulas has 
been described in numerous publications. [5, 22, 23] Many synthesis 
formulas consider the N-dimensional dynamic system of Eq. (2.1) 

[M]x(t) + [c]x(t) + [K]x(t) - f(t) (3.1) 


subject to assumptions concerning the form of the damping matrix [C], 
Perhaps the most familiar simplification stipulates that the damping 
matrix may be expressed as a linear combination of [M] and (KJ. This 
proportional damping assumption enables the classical expression for 
complex frequency response functions (w) to be written [8] 




x i (w) 

W 


N 

E 

r=l 


!r fi) !r (j) 


yi - £") + J2C r <2-)] 


(3.2) 
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The N-vectors and parameters are the normal modes and undamped 
natural frequencies calculated from the eigenvalue problem 

{ [K] - v^[M]H r = 0 (3.3) 

The scalars and are defined in Eqs. (3.4A) and (3.4B). 


K 

r 

- ^ [K] !r 

(3.4A) 


* r [c]* r 


<r 


(3.4B) 

M 

r 


(3. 4C) 


A less restrictive assertion that Is frequently encountered 
maintains that the damping matrix fC] Is symmetric. Based on this 
supposition, reference f 20 ] derives 


x (w) 2N X (J+N)X (i+N) 

V“> - 17S0 ' 4 , - T r r ~ 


(3.5) 


Equation (3.5) is based upon the state vector form of Eq. (3.1) given in 
Eqs. (2.3), (2.4), (2.5) and (2.6). The 2N vectors X r of length 2N are 
the right eigenvectors, and the 2N scalars X are the eigenvalues 
associated with Eq. (2.3), i.e.. 


U r fA] + [B]}X. - 0 


(3.6) 


The parameter a^ is the right eigenvector contraction of the matrix [A]. 


a 

r 


x t [a]x 

_r _r 


(3.7) 
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Using no simplifications regarding the form of the damping 
matrix [C], references [5] and [23] present a frequency response 
synthesis formula suitable for general linear, time- invariant systems. 


2N A. 
H . . (w) - Z 


r=l jw - A 


(3.8) 


The above representation of (w) stems from the Laplace transform of 
Eq. (3.1), 


(s [M] + s[C] + [K]} x(s) * f(s) 


(3.9A) 


x(s) - (s 2 [M] + s[C] + [K]}" l f(s) - [H(s) ] f (s) (3.9B) 


When a partial fraction expansion is carried out for each term in the 
matrix [H(s)], the form given in Eq. (3.8) results. For application in 
this thesis, however, it will be more convenient to represent FRF's for 
linear, time-invariant systems in terms of an eigenvector expansion 
similar to Eqs. (3.2) or (3.5). 


3.2 Synthesis for Linear, Time-Invariant Systems 

A frequency response synthesis formula applicable to general 
linear, time- invariant systems and expressed in terms of system eigen- 
vectors can be achieved by starting with the state vector form of 
the governing equations in Eq. (2.3). 

[A]X(t) + [B]X(t) = F(t) (3.10) 

The right and left eigenproblems associated with Eq. (3.10) are shown in 
[12] to be 
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U r fA] + rB] }X r - 0 (3. 11 A) 

U r U] T + tB] T }Y r - 0 (3.11B) 

Moreover, it is also proven that these system eigenvectors are 
bi-orthogonal with respect to the matrices [A] and [B]. 



(3.12A) 


(3.12B) 


This attribute of X r and Y f can be used to uncouple the system in Eq. 

(3.10). The solution X(t) may be expressed as a superposition of the 

<*• 

right eigenvectors. 


2N 

X(t) - Z X q (t) 
r-1 ' 


[Xjq(t) 


(3.13) 


The right eigenvectors X r comprise the columns of the 2N * 2N matrix 
[X r ] and q(t) is a 2N-vector. When Eq. (3.13) is substituted into Eq. 
(3.10), and the resulting system is premultiplied by 
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Eqs. (3.12A) and (3.12B) assure that the resulting system is 
diagonalized. 


O r ^]q(t) + j> r .Jq(t) - [Y r ] T F(t) (3.15) 

When Eq. (3.15) is Fourier transformed, a set of algebraic equations 
result. 


{jwf'a^] + ['b^] }q(w) - [Y 1 ] T F(w) 


(3.16) 


Each equation in the 2N rows of Eq. (3.16) may now be solved 
individually. The kth unknown q^(w) becomes 


<l r (w) 


Y T F(w) 

jwa + b 
J r r 


(3.17) 


The original solution X(w) is recovered when Eqs. (3.17) and (3.13) are 
combined. 


X(w) 


2N 

E 

r-1 


X Y T F(w) 
,r _r _ 

jwa + b 
J r r 


(3.18) 


Equation (3.18) assumes a more familiar form when X(w) is 
expressed in terms of the system elvenvalues A y . Premultiplication of 
the rth right eigenproblem Eq. (3.11A) by the rth left eigenvector 
generates the identity 

A Y T [A]X > . + Y^[B]X - 0 (3.19) 

r jv 


or using (3.12A) and (3.12B) 
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X a + b - 0 (3.20) 

r r r 

Substitution of Eq. (3.20) into Eq. (3.18) produces a form explicit in 
the eigenvalues X^. 

2H X Y T P(u) 

*<»> - — <3» r : x ) < 3 - 21 > 


The expansion in Eq. (3.21) provides the desired frequency 
synthesis formula, expressed in terms of eigenvectors, for any linear, 
time-invariant system. A comparison of Eqs. (3.21) and (3.8) shows that 
each FRF implied in Eq. (3.21) is given by 



(w) 


x t (w) 

fj (w) 


2N 

l 

r-1 



where the complex residue is defined in Eq. (3.23). 



X r (i+N)Y r (J+N) 


a 

r 


(3.22) 


(3.23) 


3.3 Synthesis for the Reduced System 

The parameter extraction scheme outlined in Chapter 2 provides 
approximations for the right eigenvectors X f and eigenvalues X r of the 
original system. Direct implementation of Eqs. (3.22) and (3.23) would 
additionally require the estimation of Y f and a^. However, frequency 
response function synthesis can also be accomplished by considering 
only reduced system parameters. 
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Following exactly the same steps in section (3.2) for the 
reduced system model in Eq. (2.27), 

[A] f (t) + [B]T(t) = F(t) (3.24) 

an identity similar to Eq. (3.21) may be cited for the reduced frequency 
response 

2M X Y*?( w) 

r(w) - Z - ~ r - r ~ *— (3.25) 

r°l a f (jw - X f ) 

The quantities comprising Eq. (3.25) are analogous to those comprising 
Eq. (3.21), but hold for the reduced system in Eq. (3.24) instead of Eq. 

A A A 

(3.10). X r , Y r « and X f are the right eigenvectors, left eigenvectors, 
and eigenvalues of the reduced system. 

{X r [A] + [B])X r = 0 (3.26A) 

(X r [A] T + [B] T )Y r - 0 (3.26B) 

A 

The scalar a^ is the rth left /right eigenvector condensation of the 

A 

matrix [A] . 

A AIII A A 

a = Y [A]X (3.27) 

r ,r _r 

Because all of the constituents in Eq. (3.25) can be defined from the 
parameter extraction procedure in Chapter 2, it can serve as the basis 
for generating frequency response functions. 
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Before Eq. (3.25) can be utilized to render a frequency domain 
solution to the original system of Eq. (2.3), the reconstructing 


X(w) 


in Eq. 

(2.32) 

x(w) 

r 

x(w) 



f 4> r ] [0] 
[0] [* R ] 



Y(w) 


Y(w) 


ry r(w) 0 . 28 ) 


Inserting Eq. (3.25) into Eq. (3.28), the full 2N * 1 frequency response 
X(w) is approximately reconstructed from the reduced 2M ^ 1 response 
T(w). 


A A 


2M X Y F (w) 

x( W ) S ry r • 

^ a (jw - \ ) 

r J r 


(3.29) 


When the definition of F(t) supplied in Eq. (2.27) is inserted in Eq. 
(3.29), a relationship among response and excitation spectra emerges. 


2M 

X(w) a [Y] Z 
r-1 


a Am 

X Y 
-r.r 


[D] 

[0] 


p(v) 


a (jw - \ ) 
r J r 


(3.30A) 


X(w) s [H (w) ]p (w) 


(3.30B) 


The new 2N * N matrix [H(w)] may be partitioned into N * N matrices 
P P 

* A 

[H^( w) ] and [Hjj(w)] corresponding to the velocity and displacement 
portions of the vector X(w) . 


X(w) - 

x(w) 


’nyv)] " 


x(w) 


fHpCw)] 


p(w) 


(3.31) 
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If the lower partition associated with the displacement degrees of 
freedom is extracted from Eq. (3.31), we have 

x(w) - [HpCwJJpCw) (3.32) 

It is evident that N entries in each of the columns of [Hp(w)] 
approximate the FRF's corresponding to a given exciter location. 

Although it would be more desirable to obtain all N x N 
possible FRF's, instead of only N x N^, Eq. (3.32) is sufficient for 
purposes of verification. Currently, digital Fourier analyzers natural- 
ly provide FRF's between monitored response and input locations. Thus, 
at most N x N experimental FRF's should be available for validation. 



Chapter 4 


SELECTION OF TRANSFORMATIONS 

It has been noted in section (2.1) that the accuracy achieved 
in modal parameters for a particular reduced model hinges upon the 
selection of and Two specific reduction techniques are 

considered in this chapter. Both methods attempt to automate the 
process of model reduction. 

4. 1 Independent Coordinate Method 

A widely used reduction technique common in both component 
mode synthesis [8] and experimental modal analysis T3, 7] assumes that a 
subset of coordinates may be expressed as a linear combination of the 
remaining coordinates. 


x d (t) - 


(4.1) 


In Eq. (4.1) x, (t) Is a d-vector of dependent coordinates, x. (t) is an 

-1 

i-vector of Independent coordinates, and is a d * i matrix. To 

al 

simplify notation, let x^(t) reside in the upper partition of the full, 
original response N-vector. 


x(t) 


x^t) 

x d (t) 


(4.2) 


Because it is assumed that only the independent coordinates 
need be included in the analysis, x^(t) may be immediately identified as 
y(t) in Eq. (2.9) . 
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x ± (t) = y(t) 


[# c ]x(t) 


r*c ] 


x 1 (t) 

x d ft) 


(4.3) 


As a consequence, [i^] in Eq. (4.3) may be written 


r* c i ** rm [0]] 


(4.4) 


i x N 


In Eq. (4.4), f I I is an i x i matrix and [0] is an i x d null matrix. 

The reconstructing transformation in equation Eq. (2.10) is 
also of a particularly simple form. For convenience, recall Eq. 
(2.10B). 

x(t) - [<|) R ]y(t) + e(t) (4.5) 


Using Eq. (4.1) and the fact that y(t) » x^t) produces 


Xi ( t) 

x d (t) 



f^dil 


x d (t) + e(t) 


Comparing Eqs. (4.5) and (4.6) suggests that 


(4.6) 


(*»] 



L ""di 1 J 


(4.7) 


where [I] is of i x i order. 

The only as of yet undefined quantity in [^] or |> R ] is the 

submatrix [A..J. This matrix is calculated in a least-squares sense 
al 

from measured acceleration spectra. By differentiating Eq. (4.1) twice, 
Fourier-transforming and expanding for N w discrete frequencies, it is 


possible to write 



where 


[x (w.)] «= [x, (w.) x (w ) ... x. (w )] 

' ' “ “ w i x N 


w 


and 


rx d (w k )] - [x^Wj) x d^ W 2 ^ *•* X d fw N ^ 

" “ * ‘ ~ ~ w d x N 


(4.8B) 


(4. 8C) 


When Eq. (4.8) is transposed and separated Into real and imaginary 
parts, a standard least-squares formulation emerges for [i|>,.1. 

<31 


’ m r •• i 
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[ ?i ( V J 


*di 


h ( v j 

im r .. -i 



i x d 

IM '| T 

1 - 1<Vk> J 

N 

x i 


L -J ( V J 


W V 

Upon solution of Eq. (4.9) for » Eqs. (4.4) and (4.7) completely 

determine the condensing and reconstructing transformations. 

While Blair [3] and Coppolino [ 7 J- advocate reduction methods 
employing a subset of the original coordinates, they rely on the judg- 
ment of the analyst to select the Independent and dependent coordinate 
groups. However, automatic selection procedures exist. [16] A 
procedure slightly more rigorous than manual selection is used in this 
thesis. It performs Gaussian elimination with row and column pivoting 
to select the dependent coordinate group. 

Essentially, the method permutes the rows and columns of the 


sample spectra matrix 
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[x(w i )] 


XjCwj) ... Xj(W N ) 


x,^) ... * 2 ^ W n ^ 


(4.10A) 


Xjj (w j ) ... Xjj (W^ ) 

w 


as it drives it to the form 


xxxxxxxx 

X X X X X X X 

X X X X X X 

X X X X X 

X X X X 

XXX 
X X 


ETC. 


J 



(4.10B) 


In so doing, nearly-linearly-dependent rows are driven to the bottom d 
positions shown in Eq. (4.10B). Row permutation vectors record which of 
the original rows, or variables, have been shifted to these positions. 
These variables are selected as the dependent group. This method of 
selecting the transformations will be referred to as the "independent 
coordinate method." 

When the aforementioned coordinate transformation is utilized, 
the extraction equations in section (2.3) assume a noteworthy form. 
Because the generalized, reduced coordinates are simply a subset, x^(t), 
of the original coordinates, Eq. (2.22) may be written 
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[tC] [K] fD]] 


Mx(2M+N ) 
P 


[ i£- * <v i )1 


[ — r x(w )] 

— w ^ 
i 

t-p(w i )] 


t-Y^)] 


MxN 


w 


(2M+N )xN 
p w 


(4.11) 


The reduced system equation, Eq. (2.27), and eigenvector problem, Eq. 
(2.28), may also be expressed in terms of a truncated set of the 
original response degrees of freedom. 
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*i 


F(t) 


(4.12) 


(4.13) 


Equations (4.11), (4.12) and (4.13) are precisely those utilized in the 
method of Blair and Craig. [3] 

4.2 Principal Component Method 

One distinct advantage of the method outlined in the previous 
section is that it may be interpreted intuitively. The process of 
reducing the problem size may be likened to simply picking fewer 
coordinates. Yet, many mathematically well-defined and robust reduction 
methods have been developed in the field of statistics. 

Principal component analysis is one specific reduction 
technique suggested in time series analysis literature. Fundamentally, 
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the procedure seeks to minimize in a statistical sense the error e(t) 
incurred in the reconstructed signals in Eq. (2.10B). 

Once again, development of the method starts with Eqs. (2.9) 
and (2.10). 


Y(t) - fij> c ]x(t) (4.14) 

x(t) » fi|» R ]Y(t) + e (t) (4.15) 

Principal component analysis seeks the transformations and 

such that the expectation of the inner product of the error in Eq. 
(4.15), 

E{e T (t)e(t)} (4.16) 


is minimized. For a zero mean original signal, E{x(t)} « 0, the 

transformations and ['I'n] are shown in [4] to be 

1* K 


tv - 


T 

W 

.1 


“2 


W, 


M 


(4.17) 


and 

»„> - 


(4.18) 


The vectors W comprising [0,,] and are the M eigenvectors 

«1 C K 

corresponding to the largest M eigenvalues of the variance matrix for 
x(t). The bar over [i^] in Eq. (4.18) denotes conjugation. 
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In Che analysis at hand, complex acceleration spectra are 
available. If Eqs. (4.14) and (4.15) are differentiated twice, the 
development of the principal component method outlined in Eqs. (4.16), 
(4.17) and (4.18) is unchanged. The variance matrix TS^] for the 
complex acceleration spectra may be approximated by the product 

fx(v )] [x(v )] T 

[SI s ^—1 I — i (4.19) 

N*N 

where 

[x(w A )J - [x(Wj) x(w 2 ) ... x(w N )] (4.20) 

~ ~ w 

Because [S^] is a Hermitlan matrix, it has real eigenvalues and complex 

eigenvectors. Consequently, [t|>_] and are complex when formed using 

C R 

principal component analysis. The matrices ftj> c ] and [t|» R ] are now 
completely defined in Eqs. (4.17), (4.18) and (4.19). 

Besides having a sound mathematical basis, principal component 
analysis has another attractive feature. The error measure expressed in 
Eq. (4.16) is directly related to the magnitude of the eigenvalues 
associated with the eigenvectors not retained in or [^] . 

N 

min E{e T (t)e(t)} ■ E y. (4.21) 

k-M+1 

Equation (4.21) implies that inspecting the eigenvalues of of the 
variance matrix of the acceleration spectra can indicate an appropriate 
order for the reduced model. If a distinct drop in magnitude occurs in 
the (M + l)st eigenvalue, then the first M eigenvectors of the variance 
matrix should be used in and This ability to estimate the 

number of active modes in a frequency range has no parallel in the 
independent coordinate method. 



Chapter 5 

SOLUTION OF LEAST SOUARES PROBLEMS 

Least-squares analysis is becoming more prevalent in litera- 
ture dealing vlth modal parameter identification. The methods of Blair 
[3], Coppolino [7], and Leuridan [15] are a few examples of new methods 
that use a least squares, or pseudoinverse, solution procedure. This 
thesis requires use of the method in two Instances: the solution for 

AAA 

the reduced system matrices [C], [K] , [D] and the definition of in 

al 

the independent coordinate transformation. 

The current chapter describes the essential features of the 
least-squares solution. A brief background is provided, as well as 
descriptions of two computational techniques. 

5.1 Preliminaries 

The need for a least-squares analysis can arise in the 
solution of the linear set of equations 

[Aly - b (5.1) 

The coefficient matrix [A] of the above linear system is M * N, the 
unknown vector y is N * 1, and b is M * 1. In perhaps the most common 
case, and also the problem at hand, there are more equations in Eq. 

(5.1) than there are unknowns. The overdetermined system, with M > N, 
generally need not be consistent. 

Because an exact solution might not be available, the 
classical formulation of a least-squares solution seeks to minimize the 
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L 2 norm of the residual r - b - [A]y. As will be seen, the set of 

A 

vectors, {y} that minimize the norm of the residual may contain a single 
vector, or several vectors, depending upon the specific overdetermined 
system. 

{y} ■ {y € IR N , y min I [ A] y - b |} (5.2) 

There exist other forms of least-squares analysis based on minimization 
of 1*2 and norms. [9] 

An important trait of least-squares solutions is that the 
residual r is orthogonal to the range of [A]. [14] The matrix 
formulation of this property leads to the normal equations 

[A] T (b - [A]y) - 0 (5.3A) 

or 

[A] T [A]y - [A] T b (5.3B) 

T 

When [A] is of full column rank, the coefficient matrix [A] [A] is 

positive definite and may be inverted to obtain y in Eq. (5.3B). 

A serious drawback to the normal equations method is that it 

is not as accurate for a fixed hardware word length as other methods. 

Its computational cost may be compared to other methods in terms of 

flops. A flop is the effort required to perform a floating point add, a 

multiply and some variable indexing. The normal equations method 
2 3 

requires about MN /2 + N /6 flops. [9] 

In more sophisticated approaches, orthogonal matrices play a 

A 

central role in the solution for {y} in Eq. (5.2). By definition, the 
inverse of an orthogonal matrix [Q] is its transpose. 
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[Q] T fQ] - [I] (5.4A) 

Another useful property of orthogonal matrices Is that they preserve 
Euclidean length. f21] 

||rq]x || 2 - (fQ]x) T [Q]x) « x T [0] T [Q]x - x T x = |jx|j 2 (5.4B) 

With this property, the norm in Eq. (5.2) may be replaced by an 
equivalent norm. 

II r || - ||fA]y - b|| - ||[Q] T [A]y - [0] T b|| (5.5) 

T 

A judicious choice of [Q] In Eq. (5.5) can lead to a more simplified 

expression for the norm ||r||. One way to characterize the different 

least squares techniques that follow is by the particular orthogonal 
T 

matrix [Q] they employ. 

5.2 Singular Value Decomposition 

One possible choice for [Q] in Eq. (5.5) is suggested by the 
singular value decomposition of fA], Any M x N matrix [A] may be 
expressed as the product 

[A] - [U][E][V] T (5.6) 

where 

[U] * an M x M orthogonal matrix 

[ E] ■ an M x N diagonal matrix 

[V] * an N x N orthogonal matrix 

The diagonal entries in ( E] , or singular values, are non-negative by 

definition and may be arranged in a nonincreasing order. Substituting 
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T T 

the singular value decomposition into Eq. (5.5) and choosing [Q] - [U]" 

yields 

||r|| - |imfV] T y - [U] T bj| (5.7) 


A good deal of useful information can be extracted from Eq. 
(5.7) if the following notations are introduced: 


m 




fa K J [0] 
[ 0 ] [ 0 ] 


( 5 . 8 A) 


rv] T y 


T 

[ujS 


VkJ 

§K 

?M-K 


(5.8B) 


(5.8C) 


Particular attention should be paid to the form of [E]. Only K £ N < M 
of the singular values are nonzero, which allows for N - K linearly 
dependent columns in TA]. When Eqs. (5.8A), (5.8B) and (5.8C) are 
entered into Eq. (5.7), the original norm becomes 



(5.9) 
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A straightforward application of the definition of an L ^ norm can 
simplify Eq. (5.9). If a is any vector that is partitioned into 
sub-vectors , 




' a i 


_ 1 

!i 

!i " 
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(5.10A) 


then the square of the norm of a is given by 


(a^ + ^ ... a^ ) 


(5.10B) 


II a|| ' 


(a. + at 




,2 


1 




(5.10C) 


I! a|| 1 2 - IfajH 2 + ||a 2 || 2 (5.10D) 

When the trivial property shown in Eq. (5.10D) is applied to Eq. (5.9), 
the norm of the residual may be expressed as a sum. 


||r|| 2 - ||[ a K Ihj, - gj 2 ♦ [|8 m_ k I| 2 (5.11) 


Now, one must carefully consider the rightmost term in Eq. 

(5.11). The vector g^ ^ is defined in Eq. (5.8C) to be the last M - K 

T 

entries in the product [U] b ■ g. However, [U] is uniquely defined in 

the singular value decomposition in Eq. (5.6) and depends solely upon 

T 

the matrix [A], Consequently, the product TU] b * g is fixed for a 
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given problem, i.e., a given choice of [A] and b. The selection of y 
has no affect upon flg^ ^||. 

Because Ilgj^I i® invariant with respect to any solution 

A 

vector y, the condition for minimization of the norm in Eq. (5.11) must 
be 

lit o k ]h R - gj 2 - 0 (5.12) 

From its definition, the norm of a vector is zero if and only if the 
vector itself is the zero vector 


ho^Jl^ - g R - 0 (5.13A) 

or 

C-o k .Jh K - g R (5.13B) 

Combining Eqs. (5.13B) and (5.8B) provides the least-squares solution 
sought. 


y 



(5.14A) 


• K v u^b N 

y - Z — 3 - ~ 3 ~ + E v h (5.1AB) 

J-l J-K+l ~ 3 3 


Two important conclusions can be drawn from this analysis. 

The first holds only when [A] is not of full column rank, or K < N. In 
this case, the choice of h^^ in Eq. (5.14) has no effect upon ||r|| in 

A 

Eq. (5.9). As a result, y in Eq. (5.14B) is not unique. Any choice of 
results in the same residual (fr||. In practice, h^_£ is set equal 



to 0. 
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The second conclusion pertains when [A] is of full column 
rank, or K - N. In this event, h^_^ does not exist. No ambiguity then 

A 

remains in Eq. (5.14B). The least-square solution y is unique. 

Equation (5.14B) defines the least squares solution explicitly 
in terms of the matrices obtained in a singular value decomposition of 
[A] . The actual numerical technique that produces the decomposition is 
quite lengthy. It involves application of several transformations to 
drive T A ] to bldiagonal form, then uses an iterative technique to 
achieve the final diagonal singular values. A more efficient form of 

singular value decomposition for solving least squares problems does not 

T 2 3 

form [U] and [V] explicitly. It requires approximately 2MN + 4N 

flops. [9] 


5.3 Pouseholder Transformations 

If It is known beforehand that [Al is of full column rank. 
Householder's method may be used instead of a computationally expensive 
singular value decomposition. In Householder's method, the norm 

IMI - || [Q] T [A]y - [Q] T b|| (5.15) 

is simplified by choosing [Q] to be defined from the QR decomposition 
of [A]. 


[A] - [Q][R] (5.16) 

In Eq. (5.16), [Q] is an M x M orthogonal matrix and [R] is an M x N 
upper triangular matrix. 



45 


[R] 


X X X X X X X 

X X X X X X 

X X X X X 

X X X X 

XXX 
X X 
X 


[ 0 ] 


ro] 


(5.17) 


All M * N matrices (A] may be factored as shown in Eqs. (5.16) and 
(5.17). The matrix [R^] is the M x M upper triangular part of [R]. The 
assumption that [A] is of full column rank Insures that none of the 
pivot values, or diagonal entries, of [R u ] are zero. 

Substitution of (5.16) into (5.15) and introduction of the 

notation 


[Q] T b 


§N 

. ?M-N . 


(5.18) 


generates a norm quite similar to that encountered in section (5.2). 





[R ]y 

u i 


?N 


0 


,?M-N . 



- Is 

?M-N 


(5.19) 


||r|| 2 - f fR u ly - gjj I 2 + | R m _ n II 2 (5.20) 

Following the same reasoning used in section (5.2), the condition for 
minimization of the norm in Eqs. (5.19) and (5.20) is seen to be 

[Rjy - & (5.21) 


Solution of Eq. (5.21) is particularly easy as it corresponds to the 
backsubstltutlon phase of an LU factorization. Equations (5.16), 
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(5.18) and (5.21) completely determine the solution of a least-squares 
problem via the Householder method. Because the matrix fAl Is assumed 
to be full rank, the solution In Eq. (5.21) is always unique. 

The factorization in Eq. (5.16) is never carried out 

T 

explicitly. Instead, a sequence of Householder transformations , , 

are applied to [A] to drive it to upper triangular form. A Householder 
transformation is any matrix of the form 

ro/ - fl] - 2V T V (5.22) 

T 

where V V * 1. Matrices of this type have the property that given any 

T 

two vectors x and y of equal length, a matrix [Q^]“ can be found such 
that 

tQ i ]T * " l (5.23) 

T 

The desired [Q^] is calculated using (5.22) and 

V - (x — y)/j]x — y |j (5.24) 


The diagonalization of [A] is achieved by applying a 
Householder transformation to [A] and b for each of the N columns of 
fAj. In the first transformation, x is chosen to be the first column of 


[A] and y is selected to be 



(5.25) 
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T 

When [Q^] corresponding to the above choice of x and y is applied to 
[A] and b, the first column of [A] is put in upper diagonal form. In 
the Kth subsequent step, x is chosen to be the Kth column of [A], and y 
is selected to be 




K-l.K 

o K-l _ , 

if|a K | ' >1 


(5.26) 


T 

Application of [Q ± J , calculated from the above x and y, to [Al and b 

puts the Kth column of [A] in upper triangular form, while leaving the 

previous K-l columns unchanged. The entire product of individual 

T 

transformations form the matrix [Q] in Eq. (5.15). 


([Q n 1 T [Q n _ 1 3 T ... [Qj^] 1 ) [A] - [R] 

[qi t 


(5.27) 


The computational cost of a Householder solution of the least 

squares problem is significantly less than that of singular value 

decomposition. The flop count for Householder orthogonalization is 
2 13 

MN - — N . Both singular value decomposition and Householder 
orthogonalization are of comparable stability and yield roughly twice 
as many precise digits for a given hardware word length as the normal 
equations method. [14] 



Chapter 6 

APPLICATIONS TO COMPONENT MODE SYNTHESIS 
6. 1 Preliminaries 

Often in computational dynamics, it is necessary to achieve a 
solution to the analytical model in Eq. (2.1) for a system that has many 
degrees of freedom N. 

[M]x(t) + [c]x(t) + [Kjx(t) - f(t) (6.1) 

One class of solution technique applicable to Eq. (6.1) involves inte- 
gration of the set of equations in a time-stepping procedure [2, 8]. 
However, the direct integration of Eq. (6.1) for very large N can be 
computationally expensive. 

An alternative to integration of the full system of equations 
in Eq. (6.1) is the component mode synthesis method. In addition to 
extremely large problems, this approach is also well suited to systems 
composed of a set of natural components. In the technique, the full 
system is first subdivided into substructures. Each component then 
undergoes a dynamic analysis to determine its "component modes." 

Originally, component modes were simply the free vibration 
mode shapes of the substructure, calculated with the component bound- 
aries to other components either fixed or free. In current variants of 
component mode synthesis, the component modes have come to include other 
fundamental component shapes. The deformation pattern associated with a 
unit static displacement is known as a constraint mode, while that shape 
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resulting from a unit statically applied force is known as an attachment 
mode. 

The final step in component mode synthesis involves using a 
reduced number of the component modes as Ritz vectors in a reduction 
scheme similar to Eq. (2.9). A reduced order system model is assembled 
and solved using the Ritz approximations for each component. 

In reference [24], Rubin insists that the most worthwhile 
component mode synthesis techniques should make possible the incorpora- 
tion of experimental substructure test results. The experimental 
parameter extraction scheme developed in this thesis can complement the 
component mode synthesis methods described in Howsman [12] or Chung [6], 
These particular techniques employ free-interface modes of vibration and 
a form of attachment modes as substructure Ritz vectors. 

The remainder of Chapter 6 describes how the present formula- 
tion of parameter extraction can be used to obtain experimental compo- 
nent modes to supplement, replace or simply verify analytically- 
determined free Interface mode shapes and attachment modes. Section 6.2 
briefly outlines the component mode synthesis procedure of Howsman [12], 
while section 6.3 describes the analytic formulation of the substructure 
Ritz vectors. The calculation of the appropriate Ritz vectors from 
experimental data is summarized in section 6.4. 

6.2 Component Mode Synthesis for Linear, Time-Invariant Systems 

If the system of Eq. (6.1) is comprised of two subsystems, the 
dynamics of each component may be represented in the form 
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fM a ] *a (t) + [C a ] ^a (t) + [K a ]: !a (t) " {a (t) (6 * 2A) 

[Mg]Xg(t) + [Cg]Xg(t) + fKg]Xg(t) = fg(t) (6.2B) 

where the <»-matrices are of N * N and the 0-matrices are of order 

a a 

Ng * Ng. Generally, the component mass, damping and stiffness matrices 
are obtained analytically, perhaps using finite element analysis. 
Although only two components, a and B, are described above, extension of 
the following procedure to an arbitrary number of substructures is 
straightforward. In keeping with the philosophy of accommodating the 
most general structures, Eqs. (6.2A) and (6.2B) may also be written 
using a state vector formulation 


[A a ]X«(t) ♦ [B a ]X a (t) - F 0 (t) 

[Ag]Xg(t) + [B g ]Xg(t) - Fg(t) 


(6.3A) 

(6.3B) 


X Q (t) and Xg(t) are of length 2N a and 2Ng, respectively. The state for 

mulation matrices [A 1 and [B ] are of 2N * 2N order and partitioned 

a o a a r 

as 


and 


to] r*g 

fM ] [C ] 
a J a J 


( 6 . 4 A) 


[B ] 
1 a 


-[M a ] [0] 
[0] [K a J 


C6.4B) 


The 2N_ * 2N„ matrices [A.] and [B fl ] are similarly organized. 

P P P 0 
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The components In Eqs. (6.2A) and (6.2B) must be coupled 
together to Insure that they act as a system. Conceptually, this 
condition is equivalent to requiring that the boundary response of 
component a be Identical to the adjacent boundary response of component 
6. The actual enforcement of displacement or velocity constraints is 
achieved with a locator matrix [E] that consists only of zeroes and 
ones. The matrix [E] selects and orders interface degrees of freedom in 
X(t) for each component so that 

([ElX(t)) a - ([E]X(t)) g - 0 (6.5) 

completely defines compatibility between a and 6. An additional 
constraint is that the forces exerted on substructure a by substructure 
3 must be equal and opposite to forces exerted on 6 and a. 

Fowsman [12] and Hale [11] achieve a systematic implementation 
of these constraints by using a variational principal. The details are 
quite lengthy, but a noteworthy consequence of the variational 
formulation is that the adjoint, or co-state, equations corresponding to 
Eqs. (6.3A) and (6.3B) are considered in the theory. 

■ fA a jT ta (t> + fB a ]T Ia (t) " !a (t) (6 ‘ 6A) 

-[A g ] T Y g (t) + [B g ] T Y g (t) - Fg(t) (6.6B) 

The feature of component mode synthesis that permits the 
formulation of a reduced-order system model is the introduction of P.itz 
vectors. Unlike most other methods, however, references [11] and [12] 
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suggest Rltz approximations for the standard differential equations. Eq. 
(6.3), and the adjoint differential equation, Eq. (6.6). 


X(t) 


X 

E * r (t) " f**] r T (t) 


(6.7A) 


Y(t) 


N.. 


!?YKEYK (t) 


[y r Y (t) 


(6.7B) 


The vectors ^(t) and r y (t) are x 1 and N y x 1. Equations (6.7A) and 
(6.7B) represent reductions in order when N y and N^. are chosen such that 

N x « N c ( 6 . 8 A) 

N y « N c (6.8B) 

where is the number of degrees of freedom for a component. 

The fact that not all coordinates in the collective set 

r_ e (t) ■ (r_ (t) , r_ 0 (t) } are independent becomes apparent when the Ritz 
.Jtb A a AP 

approximation, Eq. (6.7A) is substituted into Eq. (6.5). 


([EH tx ] Ex (t)) a “ (fEHy ^ (t)) 3 - 0 (6.9) 


If N D is the number of displacement and velocity constraints in Eq. 
(6.9), then of the coordinates in ^.(t) are redundant. The over- 
abundance of coordinates is remedied by requiring the user to select 

N + N- - N_ independent coordinates from r_ e (t). The set of 
(X p D JvP 

coordinates r_ c (t) is then expressible as a combination of the 



53 


independent coordinate vector rv T (t) of length N + N - N . 

% XI Cl p D 


r xs (t) 


Exi (t) 

[XD (C) 


- rc^co 


( 6 . 10 ) 


The N -vector r^C*) represents the dependent coordinates in r vc (t). 

D XD % au 

The form of the (N + N ) x (N + N - N_) matrix fC_] in Eq. (6.10) is 

ap a 6 u a 

derived to be 


[ Cx l 


m 


( 6 . 11 ) 


A completely analogous argument may be used to define a similar matrix 
[Cy] corresponding to redundant coordinates in ry S (t). 

With the transformations defined in Eqs. (6.7A), (6.7B) and 
(6.11), the variational formulation suggested in T12] yields the final 
reduced order system model. 


where 

[A s ] - 

[B s l - 


[A s ]r XI (t) ♦ [Bglr^t) - F s (t) 


(6.12) 


[C Y ] 


’f?/« o 


' f A J 0 
01 



. 0 [ {r ] B . 


0 [A s ] 


1 

QQ 

1 




(6.13A) 


[Cy] 


[ ?? ] a ° 


’(».] 0 


'v. 

1 

o 

'i 

xx> 


. 0 f V . 


1 

•>? 

tx> 

1 


[ C X ] 


(6.13B) 
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F g (t) 



r* Y i a o 


V c > 



^3 

U> 

ft 


(6. 13C) 


As in any Ritz approximation, the accuracy of the reduced-order model of 
Eq. (6.12) depends upon the selection of Ritz vectors in and [i|> ] 

«. A « A 

in Eqs. (6.13A), (6.13B) and (6.13C). 


6.3 Analytical Ritz Vectors 

Two types of Ritz vectors are used in the component mode 
synthesis method described in section (6.2). These are the so-called 
free-interface substructure modes and attachment modes. Strictly 
speaking, Eqs. (6.7A) and (6.7B) suggest that two different sets of Ritz 
vectors must be specified for each component. However, references [11] 
and ri2] advocate choosing the sets to be identical to save 
computational 
costs. 

■ 'I'y (6.14) 

Free-interface substructure modes are quite common and are 
calculated from the eigenproblem for each component o and B. The 
columns of the matrices r$_] and Tip-.] corresponding to substructure 
modes are denoted as X^ in Eq. (6.15). 

U r rA] + [B]}X r = 0 (6.15) 

Implicit in Eq. (6.15) is that the substructure boundaries are free to 
displace. If the component eigenproblems are of sufficiently low order, 
it may be feasible to calculate and separately. In this case. 
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the column vectors of [i^,] corresponding to substructure modes are 
derived from In the left elgenproblem 

U r [A] T + fB] T }Y r - 0 (6.16) 

In their most fundamental form, the attachment modes may be 
defined as the displacement shapes resulting from a statically applied 
unit force at a single degree of freedom. Reference [12] generalizes 
this definition somevhat by expressing the attachment modes in a state 
vector form. The generalized pseudostatic response of the components a 
and 6 Is of the form 


[B]X - F (6.17) 

The attachment modes in the state vector formulation are then the set of 
Ritz vectors [ip.] that satisfy 

—A 

U A ] - [B] 

[1^] Is an identity matrix defining unit forces at several displacement 
degrees of freedom. From Eq. (6.18) it Is obvious that standard 
attachment modes are simply columns of [B] *. 

When a set of Ritz vectors is composed of both substructure 
modes, Eqs. (6.15) or (6.16), and attachment modes, Eqs. (6.18), some 
redundancy results. The standard attachment modes may be expressed as a 
linear combination of the substructure modes. If K is the number of 
substructure modes kept as Ritz vectors and is the number of 
component degrees of freedom. 




(6.18) 
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[B] 


-1 


K 

Z 

r-1 



2N X Y 
Z c -r-r 

r-K+1 


-r 


(6.19) 


provides an explanation of attachment modes, or columns of [B] In 
terms of substructure modes. It Is Implicitly assumed In Eq. (6.19) 

*P 

that the product a^ ■ Y*[A]X r is normalized to unity. The rightmost 
term In Eq. (6.19) represents the contribution to the attachment modes 
of higher substructure modes that are not retained in the Rltz basis. 

A modification of standard attachment modes has been designed 
to represent the effects of the discarded higher modes. These "residual 
attachment modes" may be written 


or 




2N X Y 

* - 3 * 

r-K+1 r 


(6.20A) 





K 
Z - 
r-1 


X Y 


_r _r 



(6.20B) 


It Is important to note that the residual attachment modes may be 

calculated from (6.20B) without explicit knowledge of the deleted modes, 

X and Y for 1 - K + 1 to 2N_. 

_r ,r C 

6.4 Experimental Rltz Vectors 

Experimental analogs of the analytical free interface 
substructure modes in Eq. (6.15) have already been derived in Chapter 2. 
Right eigenvectors of order 2N for each component may be 
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approximated from the experimental reduced model eigenvectors of order 
2M C x 1. 


X 

„r 




( 6 . 21 ) 


The transformation [?_] is 2N x 2L, in this case. 

K C L 

A derivation of experimental standard attachment modes can be 
carried out by considering the pseudostatic response of the reduced 
experimental model. 


[B]r «■ F 



( 6 . 22 ) 


The elements of [B] have been completely defined in the solution of the 

A 

least squares problem in Eq. (2.24). Assuming the inverse of [Bl 

exists, the pseudostatic response of the reduced system to unit forces 

at N exciter locations is given in Eq. (6.23). 

P 


[ EiE2 


Em. 



A 

[D] 

ro] 


(6.23) 


Approximations of the standard attachment modes for the original system 
can be reconstructed from Eq. (6.23). 






'ElE2 




tf R ) (Bl" 1 


A 

[D] 

[ 0 ] 


(6.24) 


Residual attachment modes may also be estimated from the 
pseudostatic response of the reduced system. The vector T in Eq. (6.22) 
may be expressed as a superposition of the reduced right eigenvectors. 
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r 


K 

E 

r-1 


X q + 
_r r 


2M 

E C 

r-k+1 


X q 
rr 


(6.25) 


The first summation in Eq. (6.25) includes modes estimated accurately in 
the extraction procedure. These modes are analogous to the kept modes 
in Eq. (6.19). The last summation in Eq. (6.25) represents discarded, 
higher frequency modes. 

When Eq. (6.25) is substituted into Eq. (6.22) and premulti- 

A 

plied by the rth left eigenvector Y^, all terms vanish except 


or 


A A A*P A 

Y [B]X q - Y F 
i r M r r 


A A*P A 

b q - Y F 
rr _r 


(6.26A) 

(6.26B) 


With the introduction of b r - '“* r a r f torn Eq. (3.30), Eq. (6.26) can be 
recombined with Eq. (6.25). 


A A*P A A A*P A 

K X Y F 2M X Y F 

„ - r - r - » „c _r.r_ 

j + £ 

r-1 -X r a r r-K+1 


(6.27) 


As in Eq. (6.19), the left and right terms in Eq. (6.27) can be 
identified as the contribution to the reduced pseudostatic response of 
kept and deleted reduced system modes. Residual attachment modes of the 
original system relative to exciter locations can be approximated 

from Eq. (6.27) by introducing the estimate to r In Eq. (6.23) and the 
reconstructing transformation f 




A aT 1 

X Y 
r r 


K 
E 

r-1 -y r 


EDI 

0 


(6.28) 



Chapter 7 


VERIFICATION OP THE METHOD 

This chapter outlines two studies conducted to validate the 
parameter extraction method Introduced in this thesis. Two characteris- 
tics of the approach are emphasized in the verification process. 
Section (7.1) focuses upon the convergence and stability properties of 
the technique in the presence of noise when either of the transforma- 
tions discussed in Chapter 4 is utilized. In section (7.2), the 
performance of the method in resolving closely-spaced modes is compared 
to other methods. 

7. 1 Stability Example 

In the first example, the 8-DOF finite element model of the 
beam-rotor assembly shown in Figure (7.1) is considered. As apparent in 
the illustration, the model makes provision for discrete damping at all 
degrees of freedom. Moreover, the two rotors in the assembly lead to 
Coriolis damping, which manifests itself in a skew-symmetric damping 
matrix. As a result, the system modes are complex. The system mass, 
damping and stiffness matrices are listed in Tables (7.1), (7.2) and 
(7.3). 

Using the aforementioned analytic model, simulated experimen- 
tal response spectra have been generated and later analyzed by the 
parameter estimation method. By comparing the resulting approximate 
modal parameters to those obtained directly from the analytic model, the 
accuracy of the technique may be gauged. 
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Actual measurements all contain some level of noise. Inasmuch 
as It is desirable to simulate physically realizable systems, random 
noise has been added to the simulated signals in proportion to the RMS 

M •• 

response signal level. If {x^ (w^)Xj (w^) ... } represents the accelera- 
tion spectrum of the jth degree of freedom, then 


“RMS 

X j 


✓ 


Ex^ (w x ) x^ (v 2 ) 
N 


(7.1) 


is the RMS response signal level used in noise calibration. Specifi- 
cally, the noise added to the Lth frequency point of the response 
spectrum for the jth degree of freedom is calculated from 


where 


yv 

s 

w 


N 


N j fw L } 


Sxf 8 V 


x (w L ) 


JVj (w L )x j (w L } 


the complex noise added to x . (v ) 

J 

a user input noise-to-signal ratio 
a random, uniform (0,1) weighting factor 
number of frequency points in the spectrum 


(7.2) 


An extensive set of runs has been carried out for the 
analytic model. The parameters comprising the mass, stiffness and 
damping matrices in Tables (7.1), (7.2) and (7.3) have been selected 
such that two moderately damped complex modes reside in the frequency 
range from 256 to 512 rad/sec. 



DAMPED FREQUENCY 

DAMPING RATIO 

MODE 1 

309.2 (RAD/SEC) 

.08231 

MODE 2 

420.6 (RAD/SEC) 

.05444 


(7.3A) 

(7.3B) 
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In the simulation, 512 frequency points are used, with a frequency 
resolution of Aw • 0.3516 rad/sec. 

Table (7.4) depicts in matrix form the particular cases 
considered in the verification. Each entry in the matrix denotes the 
actual number of runs executed for a particular set cf conditions. 


MODEL 

ORDER 

0% 

2% 

NOISE-TO- 

4% 

■SIGNAL RATIO 
6% 

8% 

10% 

8 

10 

10 

10 

10 

10 

10 

7 

10 

10 

10 

10 

10 

10 

6 

10 

10 

10 

10 

10 

10 


TABLE (7.4) 


Each test condition is executed 10 times because random excitation has 
been simulated in all studies. As shown in the table, the noise-to- 
signal ratio is increased from 0% to 10% in 2% increments for models 
comprised of 6th to 8th order matrices. Because two reduction methods 
have been suggested in Chapter 4, the entire matrix in Table (7.4) is 
repeated for both the independent coordinate and the principal component 
transformations . 

Figures (7.2) through (7.9) summarize the error incurred in 
estimating the modal parameters for the various noise levels, model 
orders and transformation types. Each point in the figures represents 
the median of the 10 runs conducted for every case in Table (7.4). The 
noise-to-signal ratio described in Eq. (7.2) is plotted along the 
abscissa in all the graphs. The normalized error for a particular type 
of modal parameter is charted along each ordinate. Equations (7.4A), 
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(7.4B) and (7.4C) define the normalized error computed for damped 
frequencies, decay rates and the ith complex mode, respectively. 


E 


w 



E 

a 



h 


(X ± - x 1 ) T (x 1 - x ± ) 

Si 


(7.4A) 

(7.4B) 

(7.4C) 


A A A 

In Eos . (7.4A) through (7.4B), the quantities w^, a, X are the 
approximated damped frequency, decay rate and complex modes, while w^, a 
and X are exact values calculated from the analytic model. 

As suggested in previous studies [3], estimates of the damped 
frequency are the most accurate of the approximated modal parameters. 
Figures (7.2) through (7.5) depict the damped frequency normalized error 
at different noise levels for the modes at 309 rad/sec and 420 rad/sec. 
Each line on the graphs indicates the error level for a specific model 
order. For example. Figure (7.2) shows that the 6th and 7th order model 
have less error in the frequency estimates than the full order model for 
high noise levels. 

Both independent coordinate and principal component transfor- 
mations yield excellent frequency estimates. Inspection of Figures 
(7.2) through (7.5) shows that the damped frequency at 309 rad/sec has 
been estimated using an independent coordinate transformation with a 
median error less than 0.48X for all noise levels. Likewise, the 
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principal component approximation of the same mode has a median error 
less than 0.80%. Similar results hold for the mode at 420 rad/sec. 
When an independent coordinate reduction is employed, the 420 rad/sec 
damped frequency is estimated vlth a maximum median error of 0.40%. A 
median error less than 1.3" is achieved when approximating the some a: da 
using a principal component transformation. 

In the current analysis, approximations to the decay rate and 
mode shape for the 309 rad /sec mode are typical of the method's perfor- 
mance for all estimated modes. Figures (7.6) and (7.7) indicate that 
the maximum median error in decay rate is from 12% to 15% for both the 
independent coordinate and principal component reduction strategies. 
However, in both cases the error in decay rate is considerably less over 
much of the noise range considered. When independent coordinates are 
used in Figure (7.6), the decay rate actually improves in accuracy (4% 
to 6%) as the noise level is increased. All but three data points 
reside below the 7% median error level in Figure (7.7), where principal 
components have been used. 

Figures (7.8) and (7.9) depict the normalized modal vector 
error incurred in estimating the 309 rad /sec mode. In the case in which 
independent coordinates have been selected, the maximum median error is 
roughly 2.80% and occurs at the 10% noise-to-signal level. The princi- 
pal component study illustrated in Figure (7.9) has a maximum normalized 
modal error 8.0%. Again, the error in estimating the modes is much less 
over the majority of the investigated noise range. 

In reviewing the performance of the technique, both transfor- 
mation methods yield accurate parameter estimates, despite the 
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introduction of relatively high noise levels. Damped frequency 

approximations are particularly good, though the examples suggest that 

the principal component estimates typically approach the exact values 

from above. The use of an independent coordinate reduction has 

generated slightly more precise approximations than chose obtained using 

principal components. This comparison emphasizes an important attribute 

of the principal component reduction strategy; the accuracy of the 

approach depends upon the accuracy of the estimated variance matrix in 

Eq. (4.19). Moreover, the precision of the variance matrix is reliant 

upon the sample record size N in Eq. (4.20). [4] 

w 

In the test cases executed in this chapter, a record length of 
512 points has been used. Preliminary studies with even fewer frequency 
points, N w - 256, have resulted in still higher error levels for the 
principal component reduction strategy. Evidently, a tradeoff exists. 
The principal component reduction method is desirable in that it can 
provide an estimate of the number of active modes in a frequency range. 
Yet, for a small sample size, the accuracy of the method may suffer. 

7.2 Mode Resolution Example 

Section (7.1) has given an example in which the current method 
is shown to be stable in the presence of significant random noise 
levels. While such a demonstration is certainly necessary, a major 
motivation for the use of any multi-shaker method is its ability to 
extract closely-spaced modes. [13] The remainder of this section 
compares the performance of the current method to that of other 
multiple-DOF techniques in resolving closely-spaced modes. 
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Towards chis goal, an experimental modal analysis of the 
dual beams shown in Figure (7.10) has been carried out using the 
current technique, the Polyreference method, and the complex exponen- 
tial method. The beams in Figure (7.10) are designed such that the 
natural frequencies of the system occur in closely- spaced pairs. 
Uncorrelated, random excitation has been applied simultaneously at 
stations 4Z+ and 5Z+, and acceleration responses have been measured at 
locations 1Z+ to 8Z+. A frequency range of 0 to 256 Hz and a sampling 
rate of 1024 samples /second have been used in data acquisition. All 
pertinent dimensions and location labels are shown in Figure (7.10), and 
Figures (7.11A) and (7.1 IB) show photographs of the actual experimental 
setup. Tables (7.5A) and (7.5B) list the equipment and channel config- 
uration used in the experiment. A synopsis of the data acquisition 
conditions is provided in Table (7.6). 

The only difference between the data acquisition conditions 
used in the current method and that used in the Polyreference or 
complex exponential techniques is the ensemble size. As evidenced in 
Table (7.6), 10 records of time history data are collected in the runs 
using the Polyreference and complex exponential algorithms. Only one 
record is collected for use in the method Introduced in this thesis. 
The primary reason for the additional records in the Polyreference and 
complex exponential methods is that they utilize frequency response 
functions in their estimation processes. To achieve FRF's sufficiently 
smooth for parameter extraction, 10 records are used. Because the 
current method uses raw spectra, only 1 record of time history data is 
required. 
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Tables (7.7A) and (7.7B) summarize the damped frequency and 
damping level for two of the modes estimated In the analyses. Four 
runs using an independent coordinate reduction strategy are listed, as 
well as four runs using a principal component reduction scheme. 
Similarly, the results of two runs are tabulated for each of the 
Polyreference and complex exponential techniques. 

The first table gives the approximations of the lower of the 
two highly coupled modes. The damped frequencies obtained in the 
Polyreference and complex exponential techniques vary from 118. A to 
120.5 hz. When Independent coordinates have been generated, the 
damped frequency estimates range from 118.6 to 119.9 hz, in close 
agreement with the other methods. Two of the estimates achieved with 
the use of principal components are 120.2 and 120.9 hz, again in close 
correlation with the other methods. As noted in the analytical 
example problem, the damped frequencies for the principal component 
reduction strategy seem to approach the actual values from above. 

As might be expected from other studies [3], the damping 
values for all techniques are more varied. In the Polyreference and 
complex exponential algorithms, damping values range from .00255 to 
.02595. Similarly, the damping levels estimated using independent 
coordinate and principal component transformations range from .00300 
to .00535, and from .00143 to .0101, respectively. 

In Table (7.7B), approximations to the higher of the two 
closely spaced modes are listed. Correlation between the Poly- 
reference method and the current technique is again quite good. The 
complex exponential algorithm, however, yields uniformly high 
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estimates of the damped frequency and damping level. This tendency may 
be attributed In part to the relatively modest number of frames In the 
ensemble. Other runs using the complex exponential method with fewer 
frames of data have shown even higher estimates of the modal para- 
meters. 

As a final observation, both Householder's technique and 
singular value decomposition have performed comparably In terms of 
accuracy and stability in these experimental example runs. The 
singular value decomposition algorithm has required slightly more user 
Input to execute. In practical analyses, extremely small computed 
singular values may degrade accuracy and should be set equal to zero. 
Selection of the threshold value below which to zero singular values 
seems to require a significant amount of intuition. For this reason. 
Householder's method can be especially well suited for those unfamil- 
iar with least squares solution procedures. 
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DUAL BEAM DIMENSIONS 




TABLE (7.5A) 


DATA ACQUISITION EQUIPMENT 


DESCRIPTION : 

GENRAD 2515 Computer Aided Test System 1 

GENRAD 16-Channel Expansion 1 

PCB Signal Conditioner, Model 483A17 1 

LING 200 Series Vibrator 2 

LING Power Amplifier, Model TPO 25 2 

PCB Accelerometers, 303A SERIES . 8 

PCB Force Transducers,208B SERIES 2 

HP 3900 Tape Recorder 1 


TABLE (7.5B) 
CHANNEL CONFIGURATION 


CHANNEL 

NO. 

SIGNAL 

TYPE 

COOR- 

DINATE 

VOLTAGE 

RANGE 

CALIBRATION 

SCALE 

ACCELEROMETER 
SERIAL NO. 

1 

Force 

4Z+ 

-8.0, 8 

.0 

1.85185E-3 

lb/mv 

SN 3967 

2 

Force 

5Z+ 

-8.0, 8 

.0 

1.84502E-3 

lb/mv 

SN 3965 

3 

Accel. 

1Z+ 

-.125, 

.125 

0.10061 

g/mv 

SN 5915 

4 

Accel . 

2Z+ 

-.125, 

.125 

0.10204 

g/mv 

SN 4743 

5 

Accel. 

3Z+ 

-.125, 

.125 

9.31099E-2 

g/mv 

SN 4344 

6 

Accel. 

4Z+ 

-.125, 

.125 

0.10152 

g/mv 

SN 5937 

7 

Accel. 

5Z+ 

-.125, 

.125 

9.78474E-2 

g/mv 

SN 5916 

8 

Accel. 

6Z+ 

-.125, 

.125 

9.31099E-2 

g/mv 

SN 4520 

9 

Accel. 

7Z+ 

-.125, 

.125 

0.10741 

g/mv 

SN 4613 

10 

Accel. 

8Z+ 

-.125, 

.125 
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TABLE (7,6) 

DATA ACQUISITION CONDITIONS 


PROPOSED METHOD : 

Sampling Frequency (f a 1/At) 

Frequency Range 

Excitation 

References 

Responses 

Ensemble Size 

Record Size 

POLYREFERENCE/COMPLEX EXPONENTIAL : 
Sampling Frequency (f * 1/ t) 
Frequency Range ...... 

Excitation 

References 

Responses 

Ensemble Size ....... 


1024 

0-256 HZ 
Random 

2, uncorrelated 
8 
1 

2048 PTS 
1024 

0-256 HZ 

Random/Hanning Window 
2, uncorrelated 
8 

10 


Record Size 


2048 PTS 




TABLE (7.7A) 

























TABLE (7.7B) 






























Chapter 8 
CONCLUSIONS 


A modal parameter estimation procedure applicable to linear, 
time-invariant systems has been presented In this thesis. In the 
technique, multiple non-coherent input excitations may be applied to the 
structure. A consistent set of modal parameters is generated In the 
algorithm. 

The technique is formulated such that no assumptions concern- 
ing the damping matrix [C] are made. Despite relatively high noise 
levels, accurate damped frequencies, decay rates and complex modes have 
been calculated in Chapter 7 for a system with both skew-symmetric and 
non-proportional damping terms. 

As with Blair's modification of the SFD algorithm, the 
current method has been shown to accurately resolve closely spaced 
modes. An experimental comparison has shown that the parameters ob- 
tained in the current method are comparable to those achieved with the 
Polyreference and complex exponential algorithms. 

Two solution procedures noted for their stability, the singu- 
lar value decomposition technique and Householder's method, have been 
evaluated in conjunction with numerous verification runs of the method. 
In practical experimental analyses, both methods have been shown to be 
accurate and stable. Householder's method has been found more favorable 
due to its faster execution time and ease of use. 
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Since many practical modal analyses involve numerous measure- 
ment locations, two different transformations are investigated to reduce 
the size of the actual problem to be solved. Both approaches attempt to 
minimize the amount of required user interaction. The Independent 
coordinate reduction method is shown to be more accurate for a small 
sample size in several example runs. The principal component reduction 
algorithm can provide an estimate of the number of modes active in a 
given frequency range. 

Further development is needed in several areas related to this 
thesis. As an example, the procedures suggested in Chapter 6 to calcu- 
late standard and residual attachment modes from experimental data 
should undergo numerical study. Another potential application of the 
method might investigate the use of the experimental reduced systems 
matrices directly in a substructure coupling procedure. 

Modifications of the method itself may also be considered. 
Some of the statistical methods mentioned in [16] could be incorporated 
into the model reduction process. Another improvement might be to 
modify the least squares problem for the reduced system matrices in Eq. 
(2.24) such that frequency response functions are used Instead of 
spectra. In the present form, increasing the size of the experimental 
samples to achieve greater accuracy necessarily creates a larger least 
squares problem for the system matrices. If frequency response func- 
tions are used, they could be smoothed by increasing the number of 
averages, or ensemble size, used in their acquisition. The least 
squares problem for reduced system matrices could then use frequency 
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response functions smoothed by averaging. Instead of Increasing the 
number of frequency points in the raw spectra. 

Finally, additional work on the development of accurate and 
easily calculable modal confidence factors for the present method should 
be pursued. A procedure that automatically rejects the most obvious 
spurious modes would be most desirable. 
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